\documentclass[11pt]{revtex4}
\usepackage{graphicx}
\usepackage{amssymb,amsmath}
\usepackage{epstopdf}
\usepackage{bm}
\DeclareGraphicsRule{.tif}{png}{.png}{`convert #1 `dirname #1`/`basename #1 .tif`.png}


%\date{}                                           % Activate to display a given date or no date


\begin{document}
\title{Probabilities...}

\maketitle

\section{One Bubble}

\subsection{Bulk Probability}
We assume that, in the bulk, the number of appearing bubbles per units of volume and time is
\begin{equation}
	J(t) = J_0 e^{-\dfrac{W(t)}{RT}}.
\end{equation}
In a given volume $V$ without the influence of any surrounding bubble, this is equivalent
to an inhomogeneous Poisson process.
Accordingly, the probability $P$ that a cavitation has occurred within the volume $V$
at time $t$ is
\begin{equation}
	P(t) = 1 - \exp\left\lbrack - \int_0^t VJ(u)\,du\right\rbrack
	 = 1 - \exp\left\lbrack - \int_0^t \lambda(u) \, du\right\rbrack.
\end{equation}

\subsection{First Cavitation Time}
To find the time when the bubble happens, we set a uniform random variable $\xi\in\rbrack0:1\lbrack$
and wait for $t_i$ such that 
$$
	1-\xi = P(t)
$$
which is equivalent to
\begin{equation}
	\Lambda(t) = \int_0^{t} \lambda(u) \, du = - \ln \xi.
\end{equation}

The inversion of $\Lambda$ to find $t$ depends on the shape of $\lambda$.



\section{A set of coupled bubbles}
\subsection{Brute Force}
Let us assume that we initially have $N_0$ slots with $\lambda_1,\ldots,\lambda_{N_0}$ intensities
starting from $t_0$. Each slot has an exponential time $\xi_i$ to reach, and
we should solve $N_0$ integrals and keep the lowest time...and discard the $N_0-1$ other results...

\subsection{Finding a cavitation time and slot}
 To find the first cavitation, we form
\begin{equation}
	\hat\lambda_0(t) = \sum_{i=1}^{N_0} \lambda_i(t)
\end{equation}
and we look for 
\begin{equation}
	\hat\Lambda_0(t) = \int_{t_0}^{t} \hat\lambda_0(u) \, du = - \ln \xi_0,
\end{equation}
which provides $t_1$.

To find which slot is activated, we need another uniform random variable $\zeta\in\rbrack0:1\lbrack$
and we choose a slot $i_0$ proportionally to their intensities. 

The method has been tested by Ogata and Lewis (1981) and exactly proved by Reynaud-Bouret and Chevalier (2013).



\subsection{Next One}
Then we remove $i_0$ from the list
and we go on.
After a renumbering of the $N_1=N_0-1$ slots, we are left to solve
$$
	\hat\Lambda_1(t) = \int_{t_1}^{t}  \hat\lambda_1(u) \, du = \int_{t_1}^{t}  \sum_{i=1}^{N_1} \lambda_i(t) \, du = - \ln \xi_1
$$

\subsection{Special case: majored intensities}
We don't need to perform the integral inversion, but we use a rejection method.
Let us assume that 
$$
	\forall t\geq t_0, \hat\lambda(t) \leq \bar\lambda_0.
$$
The we choose and exponential time $\hat t_1 = \hat t_0 -({\ln \xi_0})/{\bar\lambda}$.
The we choose another uniform random number $\zeta_0$ and locate $\bar\lambda\zeta_0$ in the sum of the individual intensities. 
If $\bar\lambda\zeta_0\geq\hat\lambda(t_1)$, then we reject the step.
We compute $\bar\lambda_1=\hat\lambda(t_1)$ and we start over...


\section{Numerical Inversion}
Let us solve the generic case
\begin{equation}
	\Lambda(t) = \int_{t_0}^{t} \lambda(u) \, du = - \ln \xi = \delta.
\end{equation}
Starting from $\tau_n$ an approximate solution, we get the Newton's sequence
$$
	\tau_{n+1} = \tau_n + \dfrac{\delta-\Lambda(\tau_n)}{\lambda(\tau_n)}
$$


\section{The One Next Bubble Problem in Avalanches}
\subsection{Description}
According to Brennen, the energy cost for one bubble is
$$
		W(t) = \frac{16\pi \gamma^3}{3(\Delta P )^2}
$$
and
$$
	J(t) = J_0 e^{-\frac{W(t)}{kT}}, \;\;J_0 = N\left(\dfrac{2\gamma}{\pi m}\right)^{1/2}.
$$
The nucleation time of one cavity with volume $V$ is given by the inversion of
$$
	\int_0^{t} \lambda_0 e^{-\dfrac{a}{(\Delta P)^2}}\,du = -\ln \xi
$$
with
$$
	\lambda_0 = V J_0 = V N\left(\dfrac{2\gamma}{\pi m}\right)^{1/2} \mbox{ in } \mathrm{s}^{-1}
$$
and 
$$
	a = \dfrac{16\pi\gamma^3}{3kT} \mbox{ in } \mathrm{Pa}^{2} \simeq (1.2\,\mathrm{GPa})^2.
$$

One a bubble appears in one cavity, let's say at $t=0$, the the pressure in the cavity is
$$
	P_A(t) = P_V + (P^\ast-P_V)\cos(\omega t)e^{-t/\tau}
$$
Where $P^\ast$ is the (negative) pressure  of nucleation.

We want to solve $y=P_B-P_V$
using the \textbf{real part} of
$$
	\ddot y + \nu \dot y + \omega_0^2 y = \delta P e^{(\j \omega-k) t}
$$
with
$$
\begin{cases}
	y       = Z e^{(\j \omega-k) t}\\
	\dot y  =  (\j \omega-k) Z e^{(\j \omega-k) t} \\ %+ \dot Z e^{(\j \omega-k) t} +
	\ddot y = (\j \omega-k)^2 Z e^{(\j \omega-k) t}
\end{cases}
$$
we have a $Z$ equation
$$
		(\j \omega-k)^2 Z + \nu (\j \omega-k) Z + \omega_0^2 Z = \delta P
$$
meaning 
$$
	Z = \dfrac{\delta P}{\omega_0^2+\nu(\j \omega -k)+ (\j \omega -k)^2}
$$
or
$$
	Z = \dfrac{\delta P}{\left(\omega_0^2-\omega^2 + k^2 -\nu k \right) + \j\omega \left(\nu-2 k\right)}
$$

\subsection{Simplified}
We assume $P(t)= P_V + (P^\ast-P_V)\cos(\omega t-\phi)e^{-t/\tau}$, $P^\ast\ll -P_V$.
We have
$$
	\Lambda(t) = \int_0^{t} \lambda_0 \exp\left\lbrack \dfrac{-\zeta}{\left\lbrack \cos(\omega u-\phi)e^{-ku}\right\rbrack^2}\right\rbrack \, du = -\ln \xi
$$
with
$$
	\zeta = \dfrac{16\pi\gamma^3}{3kT (P^*-P_V)^2}
$$
The only active values are set around the $u_n$ such that the denominator is maximum.
The maxima of
$$
	f(t) = -\dfrac{1}{\cos(\omega t-\phi)^2e^{-2kt}}
$$
are defined for
$$
	\left( \omega \sin(\omega t^\star-\phi) + k \cos(\omega t^\star-\phi) \right) \cos(\omega t^\star-\phi ) = 0
$$

%The zeros of $f$ and of $\lambda$ are defined by 
%$$
%	u_{0,n} = (2n+1)\dfrac{\pi}{2}
%$$
%The maxima of $f$ and of $\lambda$ are defined by
%$$
%	\omega \sin(\omega u^\star-\phi) + k \cos(\omega u^\star-\phi) = 0.
%$$
We are left with a sequence of spikes located at
$$
	\omega t^\star_n - \phi = n\pi - 2 \arctan\left(\dfrac{k}{\omega+\sqrt{\omega^2+k^2}}\right)
$$
and bracketed by 
$$
	\omega t^\ominus_n - \phi = n\pi - \dfrac{\pi}{2} \mbox{ and } \omega t^\oplus_n - \phi = n\pi + \dfrac{\pi}{2}
$$

%$$
%	\dfrac{1}{f^2(u)} \sim \dfrac{1}{f^2(u^\star)} - \dfrac{f''(u^\star)}{f^3(u^\star)} \left(u-u^\star\right)^2
%$$


It isn't possible to compute an accurate expression for this integral, so we are going to use the method of the steepest descent. \\ 
In other word, we have to find the maxima of the function. Then when the maximum is found. 
$$
\int_{0}^{t} \lambda_{0} \exp{\left\lbrack-\dfrac{\zeta}{\left(\cos{\left(\omega u -\phi\right)}e^{-ku}\right)^2}\right\rbrack \, }du 
$$

So we could use the method of the steepest descent. To do that, we need an interval, in which the function have a maximum. We are going to do a Taylor expansion around this maxima. 
Then to use this method, we have to divide the interval $[0,t]$ into N sub-interval in which the function has a single maximum. \\ 
In other word the integral is : 
\begin{align}
\nonumber
\int_{0}^{t} \lambda_{0} \exp{\left\lbrack-\dfrac{\zeta}{\left(\cos{\left(\omega u -\phi\right)}e^{-ku}\right)^2}\right\rbrack \, }du = \int_{0}^{t^{\ominus}_{1}}\lambda_{0}\exp{\left[-\dfrac{\zeta}{\left(\cos{\left(\omega u-\phi\right)}e^{-ku}\right)^2}\right]}+\\
\nonumber
+\sum_{n=1}^{N}\lambda_{0} \int_{t^\ominus_n}^{t^\oplus_n} \exp{\left[-\dfrac{\zeta}{\left(\cos{\left(\omega u -\phi\right)}e^{-k u}\right)^2}\right]}du +\int_{t^{\oplus}_{N}}^{t}\lambda_{0}\exp{\left[-\dfrac{\zeta}{\left(\cos{\left(\omega u-\phi\right)}e^{-ku}\right)^2}\right]}
\end{align}

with $t^{\oplus}_n$ and $t^{\ominus}_n$ are both the minima of the function on this interval. \\

Let $g(t)$ be equal to  $-\dfrac{1}{\left(\cos{\left(\omega t -\phi\right)}e^{-k t}\right)^2}$. \\
The integrant becomes $\exp{\left[\zeta g(t)\right]}$, and the integral : 
$$
\int_{0}^{t} \exp{\left[\zeta g(u)\right]}du = \sum_{n=0}^{N} \int_{t^{\ominus}_{n}}^{t^{\oplus}_{n}} e^{\zeta g(u)} du
$$
In addition, we can do a change of variable. Then we choose X such that: $X= \omega u-\phi$. We introduce a new function $f(X)$ such that $f(X) = g\left(\dfrac{X}{\omega}+\phi\right)$. 
The maxima of $f(t)$ is $t^{\star}_{n}$, then the taylor expansion around this maxima is: 
$$
\exp{\left[\zeta f(t)\right]} \approx \exp{\left[\zeta\left(f(t^{\star}_{n})+\dfrac{\left(t-t^{\star}_{n}\right)^2}{2}\partial_{t}^2 f(t^{\star}_{n})\right)\right]}
$$
The integral becomes : 
\begin{align}
\nonumber
\int_{0}^{t} e^{\zeta g(u)} du =
\dfrac{\lambda_{0}}{\omega}e^{\zeta f(X^{\star}_{0})}\left\{ \int_{-\phi}^{X^{\ominus}_{1}}\exp{\left[\dfrac{\zeta}{2}\left(X-X^{\star}_{0}\right)^2\partial^2_{X} f(X^{\star}_{0}) \right] }dX+ \int_{X^{\oplus}_{N}}^{X_{t}}\exp{\left[\dfrac{\zeta}{2}\left(X-X^{\star}_{N+1}\right)^2\partial^2_{X} f |_{X=X^{\star}_{N+1}} \right] }dX\right\}+\\ 
\nonumber
\sum_{n=1}^{N} \dfrac{\lambda_{0}}{\omega} \int_{X^{\ominus}_{n}}^{X^{\oplus}_{n}}e^{\zeta f(X^{\star}_{n})}\exp{\left[\dfrac{\zeta}{2}\left(X-X^{\star}_{n}\right)^2\partial^2_{X} f(X^{\star}_{n}) \right] }dX 
\end{align}

\begin{align}
\nonumber
\int_{0}^{t} e^{\zeta g(u)} du =
\dfrac{\lambda_{0}}{\omega}e^{\zeta f(X^{\star}_{0})}\left\{ \int_{-\phi}^{X^{\ominus}_{1}}\exp{\left[\dfrac{\zeta}{2}\left(X-X^{\star}_{n}\right)^2\partial^2_{X} f(X^{\star}_{n}) \right] }dX+ \int_{X^{\oplus}_{N}}^{X_{t}}\exp{\left[\dfrac{\zeta}{2}\left(X-X^{\star}_{n}\right)^2\partial^2_{X} f(X^{\star}_{n}) \right] }dX\right\}+\\ 
\nonumber
\sum_{n=1}^{N} \dfrac{\lambda_{0}}{\omega} \int_{-\infty}^{+\infty}e^{\zeta f(X^{\star}_{n})}\exp{\left[\dfrac{\zeta}{2}\left(X-X^{\star}_{n}\right)^2\partial^2_{X} f(X^{\star}_{n}) \right] }dX 
\end{align}


We will compute the two first integrals with mathematica, and we are going to compare with their accurate values. So we put in a tabular all values of the integrals. Let $I_{11}$ and $I_{12}$ be the exact  value of the integral for the first and second  spike respectively. Also let $I_{21}$ and  $I_{22}$ be the approximate value of the integral for the first and second  spike respectively.\\
\\
\begin{tabular}{|c | c | c | c| c | c |  c | c | c | c |}
\hline
$\zeta$ & $\phi$ & $\omega$ & $k$ & $I_{11}$ & $I_{12}$ & $\dfrac{I_{12}}{I_{11}}$ & $I_{21}$ & $I_{22}$& $\dfrac{I_{22}}{I_{21}}$\\
\hline
 1  & 0 & 1 & 0.1& 0.5 & 0.17 & 0.34 & 0.64 & 0.2 & 0.32\\
 \hline
 100  & 0 & 1 & 0.1& $1.77*10^{-44}$ & $3.26*10^{-82}$ & $1.84*10^{-38}$ & $1.78*10^{-44}$ & $3.27*10^{-82}$ & $1.83*10^{-38}$\\
 \hline
\end{tabular}
\end{document}